/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		shape_interface.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-01-17 16:54
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:42
*
*   Description:
*
=============================================================================*/
#ifndef SHAPE_INTERFACE_H
#define SHAPE_INTERFACE_H

#include "../../arithmetic/sorting.hpp"
#include "../../model/dof.hpp"
#include "../../model/gauss_point.hpp"
#include "../../util/factory.hpp"
#include "../../util/iterator.hpp"
#include "../basic_geometry.hpp"
#include "dimension_type.hpp"
#include <boost/geometry.hpp>
#include <boost/mpl/int.hpp>
#include <boost/property_tree/ptree.hpp>
#include <unordered_set>

using mpl1 = boost::mpl::int_<1>;
using mpl2 = boost::mpl::int_<2>;
using mpl3 = boost::mpl::int_<3>;

namespace geoxfem {
namespace geometry {

struct get_vertex_gauss;

template <typename DimensionType, typename degree>
struct Vertex;
template <typename DimensionType, typename degree>
struct Segment;
template <typename DimensionType, typename degree>
struct PolygonInterface;

enum class topology { /*{{{*/
                      /// no shared part
                      seprate,
                      /// share some part but no shared facet
                      touch,
                      /// share a facet. only between two item with same
                      /// topology dimension
                      adjacent,
                      /// contains
                      parent,
                      /// is part of
                      child,
                      /// same item
                      identical,
                      /// error
                      error = 999
}; /*}}}*/

template <typename DimensionType, typename degree>
struct ShapeInterface
    : std::vector<typename traits::point_type<DimensionType>::type>,
      util::Indexed {

  using Degree     = degree;
  using point_type = typename traits::point_type<DimensionType>::type;
  using PointList  = std::vector<point_type>;

  using Shape        = ShapeInterface<DimensionType, degree>;
  using ShapeList    = std::vector<Shape*>;
  using SubShape     = ShapeInterface<DimensionType, mpl1>;
  using SubShapeList = std::vector<SubShape*>;

  using factory  = dp14::factory<Shape, const std::vector<point_type>&>;
  using Material = material::MaterialInterface<DimensionType>;

  using Stress = arithmetic::Stress<DimensionType>;
  using Strain = arithmetic::Strain<DimensionType>;

  using Gauss     = model::GaussPoint<DimensionType>;
  using GaussList = std::vector<Gauss>;

  static constexpr size_t   dimension = traits::dimension<DimensionType>::value;
  static const inline char* poly_degree[] = {"\0", "Linear", "Quadratic",
                                             "Cubic"};

  /*constructors and destructor{{{*/

  /**
   * @brief normal constructor
   *
   * @param listPoint
   */
  ShapeInterface(const std::vector<point_type>& listPoint)
      : std::vector<point_type>{listPoint},
        util::Indexed{},
        m_boundary{false},
        m_list_gauss{},
        m_material{},
        m_xfem{false},
        m_parent(nullptr),
        m_list_child{},
        m_listconnectivity{4} {
    m_hashkey.clear();
    std::vector<int> list;
    for(auto& p : listPoint)
      list.emplace_back(p.getIndex());
    std::sort(list.begin(), list.end());
    for(auto& p : list)
      m_hashkey.append(std::to_string(p) + ",");
    m_hashkey.pop_back();
  }
  ShapeInterface() = default;
  virtual ~ShapeInterface() {
    for(auto v : m_list_child)
      delete v;
  }

  static Shape::factory& getFactory() {
    return Shape::factory::instance();
  }
  /*}}}*/
  /*miscs{{{*/

  /**
   * @brief override operator==
   *
   * @param rhs
   *
   * @return
   */
  bool operator==(const Shape& rhs) const {
    return m_hashkey == rhs.m_hashkey;
  }
  /**
   * @brief get Class Name
   *
   * @return
   */
  virtual std::string getClassName() const {
    return "ShapeInterface";
  }

  virtual Vertex<DimensionType, degree>& toVertex() {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  virtual const Vertex<DimensionType, degree>& toVertex() const {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  virtual Segment<DimensionType, degree>& toSegment() {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  virtual const Segment<DimensionType, degree>& toSegment() const {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  virtual PolygonInterface<DimensionType, degree>& toPolygon() {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  virtual const PolygonInterface<DimensionType, degree>& toPolygon() const {
    GEOXFEM_NOT_IMPLEMENTED();
  }
  /*}}}*/
  /*direct member access{{{*/

  /**
   * @brief getMaterial
   *
   * @return
   */
  const Material* getMaterial() const {
    return m_material;
  }
  /**
   * @brief setMaterial
   *
   * @param mat
   */
  void setMaterial(const Material* mat) {
    m_material = mat;
  }
  /**
   * @brief get dimension
   *
   * @return
   */
  static size_t getDimension() {
    return dimension;
  }
  /**
   * @brief get shape's hash key
   *
   * @return
   */
  GEOXFEM_Hash getHashKey() const {
    return m_hashkey;
  }
  /**
   * @brief getCoords
   *
   * @return
   */
  MatrixXR getCoords() const {
    size_t   nn = getNumberOfPoint();
    MatrixXR coord(dimension, nn);
    for(size_t i = 0; i != nn; ++i)
      coord.col(i) = getPoint(i).getCoord();
    return coord;
  }
  /*}}}*/
  /*iterator access{{{*/

  typename PointList::iterator getPointBegin() {
    return this->begin();
  }
  typename PointList::iterator getPointEnd() {
    return this->end();
  }
  typename PointList::const_iterator getPointBegin() const {
    return this->begin();
  }
  typename PointList::const_iterator getPointEnd() const {
    return this->end();
  }

  typename GaussList::iterator getGaussPointBegin() {
    return m_list_gauss.begin();
  }
  typename GaussList::iterator getGaussPointEnd() {
    return m_list_gauss.end();
  }
  typename GaussList::const_iterator getGaussPointBegin() const {
    return m_list_gauss.begin();
  }
  typename GaussList::const_iterator getGaussPointEnd() const {
    return m_list_gauss.end();
  }

  typename ShapeList::iterator getNeighbourBegin(const size_t& dim) {
    return m_listconnectivity.at(dim).begin();
  }
  typename ShapeList::iterator getNeighbourEnd(const size_t& dim) {
    return m_listconnectivity.at(dim).end();
  }
  typename ShapeList::const_iterator
  getNeighbourBegin(const size_t& dim) const {
    return m_listconnectivity.at(dim).begin();
  }
  typename ShapeList::const_iterator getNeighbourEnd(const size_t& dim) const {
    return m_listconnectivity.at(dim).end();
  }

  typename SubShapeList::iterator getChildBegin() {
    return m_list_child.begin();
  }
  typename SubShapeList::iterator getChildEnd() {
    return m_list_child.end();
  }
  typename SubShapeList::const_iterator getChildBegin() const {
    return m_list_child.begin();
  }
  typename SubShapeList::const_iterator getChildEnd() const {
    return m_list_child.end();
  }
  /*}}}*/
  /*point manuvar{{{*/
  /**
   * @brief getNumberofPoint
   *
   * @return
   */
  virtual size_t getNumberOfPoint() const {
    GEOXFEM_NOT_IMPLEMENTED();
    return 0;
  }
  /**
   * @brief getOrder
   *
   * @return
   */
  /* size_t getOrder() const { return order::value; } */
  /**
   * @brief &getHeadPoint
   *
   * @return
   */
  point_type& getHeadPoint() {
    return this->front();
  }
  /**
   * @brief &getHeadPoint
   *
   * @return
   */
  const point_type& getHeadPoint() const {
    return this->front();
  }
  /**
   * @brief get last of point
   *
   * @return
   */
  point_type& getTailPoint() {
    return this->back();
  }
  /**
   * @brief &getTailPoint
   *
   * @return
   */
  const point_type& getTailPoint() const {
    return this->back();
  }
  /**
   * @brief getPoint
   *
   * @param n
   *
   * @return
   */
  point_type& getPoint(const size_t& n) {
    return this->at(n);
  }
  /**
   * @brief getPoint
   *
   * @param n
   *
   * @return
   */
  const point_type& getPoint(const size_t& n) const {
    return this->at(n);
  }
  /**
   * @brief getPointIndices
   *
   * @return
   */
  std::vector<size_t> getPointIndices() const {
    std::vector<size_t> v;
    for(auto& p : *this)
      v.emplace_back(p.getIndex());
    return v;
  }

  /*}}}*/
  /*gauss point manuvar{{{*/

  /**
   * @brief buildStdGaussPoints
   */
  virtual void getStdGaussPoints(std::vector<Gauss>&) {
    return;
  }
  /**
   * @brief clear Gauss Point list
   */
  void clearGaussPoints() {
    m_list_gauss.clear();
    return;
  }
  void append(Gauss& gp, const _real& tol = static_cast<_real>(1e-7)) {
    if(gp.getWeight() > tol) {
      gp.setLocalCoord(this->local2global(gp.getCoord()));
      gp.setMaterial(m_material);
      gp.setVisibilityLength(getNumberOfNeighbors(0));
      m_list_gauss.emplace_back(std::move(gp));
    }
    return;
  }
  void append(Gauss&& gp, const _real& tol = static_cast<_real>(1e-7)) {
    if(gp.getWeight() > tol) {
      gp.setLocalCoord(this->local2global(gp.getCoords()));
      gp.setMaterial(m_material);
      gp.setVisibilityLength(getNumberOfNeighbors(0));
      m_list_gauss.emplace_back(gp);
    }
    return;
  }
  const Gauss& getGaussPoint(const size_t& n) const {
    return m_list_gauss.at(n);
  }
  void setNodalWeight() const {
    for(auto& gp : m_list_gauss) {
      VectorXR N   = getStdNByGlobal(gp.getLocalCoord());
      size_t   loc = 0;
      for(auto& vertex : neighborsIn(*this, 0))
        vertex->toVertex().getGauss().addWeight(N[loc++] * gp.getWeight());
    }
    return;
  }
  size_t getNumberOfGaussPoints() const {
    return m_list_gauss.size();
  }
  /*}}}*/
  /* topology related {{{*/

  /**
   * @brief get Topology Dimension of current shape
   *
   * @return
   */
  virtual size_t getTopologyDimension() const {
    return 0;
  }
  size_t getMeshTopologyDimension() const {
    for(size_t topo_dim = 1; topo_dim != 4; ++topo_dim)
      if(getNumberOfNeighbors(topo_dim) == 0)
        return topo_dim - 1;
    return 3;
  }
  bool ifVertex() const {
    return getTopologyDimension() == 0;
  }
  bool ifSegment() const {
    return getTopologyDimension() == 1;
  }
  bool ifPolygon() const {
    return getTopologyDimension() == 2;
  }
  bool ifCubic() const {
    return getTopologyDimension() == 3;
  }
  bool ifCell() const {
    if(ifCubic())
      return true;
    return getNumberOfNeighbors(getTopologyDimension() + 1) == 0;
  }
  bool ifFacet() const {
    return getMeshTopologyDimension() == getTopologyDimension() + 1;
  }
  void updateBoundaryStatus() {
    if(!ifFacet())
      return;
    m_boundary = (getNumberOfNeighbors(getTopologyDimension() + 1) == 1);
  }
  bool ifBoundary() const {
    if(ifFacet())
      return m_boundary;
    else
      for(auto& _facet : neighborsIn(*this, getMeshTopologyDimension() - 1))
        if(_facet->ifBoundary())
          return true;
    return false;
  }
  _real getArea() const {
    _real area = 0;
    for(auto& gp : gpsIn(*this))
      area += gp.getWeight();
    return area;
  }
  size_t getNumberOfNeighbors(const size_t& dim) const {
    return m_listconnectivity.at(dim).size();
  }
  bool hasNeighbour(const Shape& shp) const {
    auto& list = m_listconnectivity.at(shp.getTopologyDimension());
    return std::find(list.begin(), list.end(), &shp) != list.end();
  }
  void link(Shape& shp) {
    auto& list = m_listconnectivity.at(shp.getTopologyDimension());
    if(!hasNeighbour(shp))
      list.emplace_back(&shp);
  }
  void link(Shape* shp) {
    auto& list = m_listconnectivity.at(shp->getTopologyDimension());
    if(!hasNeighbour(*shp))
      list.emplace_back(shp);
  }
  void detachFrom(const Shape& shp) {
    auto& list = m_listconnectivity.at(shp.getTopologyDimension());
    auto  it   = std::find(list.begin(), list.end(), &shp);
    if(it != list.end())
      list.erase(it);
  }
  void detachFrom(const Shape* shp) {
    detachFrom(*shp);
  }
  virtual std::vector<std::unique_ptr<Shape>> createFacet() const {
    GEOXFEM_NOT_IMPLEMENTED();
    return std::vector<std::unique_ptr<Shape>>{};
  }

  /**
   * @brief getTopologyWith
   *
   * @param rhs
   *
   * @return
   */
  topology getTopologyWith(const Shape& rhs) const {
    /// shape1 is shape2
    if(*this == rhs)
      return topology::identical;

    size_t topo_dim1 = getTopologyDimension();
    size_t topo_dim2 = rhs.getTopologyDimension();
    /// get shared point
    std::vector<size_t> list1 = getPointIndices();
    std::vector<size_t> list2 = rhs.getPointIndices();
    std::sort(list1.begin(), list1.end());
    std::sort(list2.begin(), list2.end());
    std::vector<std::size_t> intersection_list;
    std::set_intersection(
      list1.begin(), list1.end(), list2.begin(), list2.end(),
      std::inserter(intersection_list, intersection_list.begin()));
    const size_t nShare = intersection_list.size();
    /// no contact between two entities
    if(nShare == 0)
      return topology::seprate;
    /// current shape is a part of shape2
    if(nShare == list1.size())
      return topology::child;
    /// current shape contains shape2
    if(nShare == list2.size())
      return topology::parent;
    /// if current shape and shape2 shared the same topology _dimension
    if(topo_dim1 == topo_dim2 && nShare > degree::value * (topo_dim1 - 1))
      return topology::adjacent;
    /// if others
    return topology::touch;
  }
  /**
   * @brief ifContains
   *
   * @param _pShape
   *
   * @return
   */
  bool ifContains(const Shape& shape) const {
    return getTopologyWith(shape) == topology::parent;
  }
  /**
   * @brief getNeighbourSiblingAt
   *
   * @param p_facet
   *
   * @return
   */
  Shape* getNeighbourSiblingSharing(const Shape& share_shp) const {
    const std::size_t topo_dim = getTopologyDimension();
    /* exception: find neighbour vertex sharing the same segment */
    auto it = getNeighbourBegin(0);
    if(topo_dim == 0) {
      GEOXFEM_ASSERT(getTopologyWith(share_shp) == topology::child);
      it = std::find_if(getNeighbourBegin(0), getNeighbourEnd(0),
                        [share_shp](Shape* p_sibling) {
                          return share_shp.ifContains(*p_sibling);
                        });
    }
    else {
      /* find neighbour sibling sharing the same facet */
      GEOXFEM_ASSERT(getTopologyWith(share_shp) == topology::parent);
      it = std::find_if(getNeighbourBegin(topo_dim), getNeighbourEnd(topo_dim),
                        [share_shp](Shape* p_sibling) {
                          return p_sibling->ifContains(share_shp);
                        });
    }
    if(it == getNeighbourEnd(topo_dim))
      return nullptr;
    else
      return *it;
  }
  /*}}}*/
  /*shape function related{{{*/

  /** calculate standard shape function using local coordinates
   *
   * @param lcoord
   * @return
   */
  virtual VectorXR getStdNByLocal(const VectorXR&) const {
    GEOXFEM_NOT_IMPLEMENTED();
    return Vector3R::Zero();
  }
  /** calculate standard shape function derivation dNdlx using local
   * coordinates
   *
   * @param lcoord
   * @return
   */
  virtual MatrixXR getStddNdlxByLocal(const VectorXR&) const {
    GEOXFEM_NOT_IMPLEMENTED();
    return Eigen::Matrix<_real, dimension, dimension>::Zero();
  }
  /** calculate standard shape function derivation dNdlx using local
   * coordinates
   *
   * @param lcoord
   * @return
   */
  MatrixXR getdNdlxByLocal(const VectorXR&) const {
    GEOXFEM_NOT_IMPLEMENTED();
    return Eigen::Matrix<_real, dimension, dimension>::Zero();
  }

  /** also you should get shape function N with global coordinate gcoord
   *
   * @param gcoord
   * @return
   */
  virtual VectorXR getStdNByGlobal(const VectorXR&) const {
    GEOXFEM_NOT_IMPLEMENTED();
    return Vector3R::Zero();
  }
  /** then global coordinate can be yield by function above
   *
   * @param gcoord
   * @return
   */
  virtual VectorXR global2local(const VectorXR&) const {
    GEOXFEM_NOT_IMPLEMENTED();
    return Eigen::Matrix<_real, dimension, 1>::Zero();
  }
  /** calculate jacobi matrix using local coord
   *
   * @param lcoord
   * @return
   */
  MatrixXR getJacobByLocal(const VectorXR& lcoord) const {
    /// jacobi = [dNdx] [coord]^T
    return getStddNdlxByLocal(lcoord) * getCoords().transpose();
  }
  /** then using dNdlx and jacobi matrix you can get shape function derivation
   * at global coordinates dNdgx
   *
   * @param lcoord
   * @return
   */
  MatrixXR getdNdgxByLocal(const VectorXR& lcoord) const {
    return getJacobByLocal(lcoord).inverse() * getStddNdlxByLocal(lcoord);
  }
  /** using standard shape function to get global coord using local coord
   *
   * @param lcoord
   * @return
   */
  VectorXR local2global(const VectorXR& lcoord) const {
    return getCoords() * getStdNByLocal(lcoord);
  }
  /** calculate standard shape function derivation dNdlx using global
   * coordinates
   *
   * @param gcoord
   * @return
   */
  MatrixXR getdNdlxByGlobal(const VectorXR& gcoord) const {
    return getdNdlxByLocal(global2local(gcoord));
  }
  /** calculate jacobi matrix using global coord
   *
   * @param gcoord
   * @return
   */
  MatrixXR getJacobByGlobal(const VectorXR& gcoord) const {
    return getJacobByLocal(global2local(gcoord));
  }
  /** then using dNdlx and jacobi matrix you can get shape function derivation
   * at global coordinates dNdgx
   *
   * @param gcoord
   * @return
   */
  MatrixXR getdNdgxByGlobal(const VectorXR& gcoord) const {
    return getdNdgxByLocal(global2local(gcoord));
  }

  /*}}}*/
  /*xfem related{{{*/

  /**
   * @brief setXfem
   *
   * @param value
   */
  virtual void setXfem(const bool& value) {
    m_xfem = value;
  }
  /**
   * @brief ifXfem
   *
   * @return
   */
  bool ifXfem() const {
    return m_xfem;
  }
  void updateGaussXfem(const Segment<DimensionType, degree>& crack_segment) {
    point_type p0 = crack_segment.getHeadPoint();
    point_type p1 = crack_segment.getTailPoint();
    for(auto& gp : gpsIn(*this)) {
      int    side = if_at_left_wing<point_type>(p0, p1, gp);
      size_t id   = 0;
      for(auto vertex : neighborsIn(*this, 0)) {
        std::vector<point_type> gp2vertex = {vertex->getHeadPoint(), gp};
        bool visible = boost::geometry::intersects(gp2vertex, crack_segment);
        if(!visible)
          gp.setInvisible(id, 2 * side);
        id++;
      }
    }
    return;
  }
  /*}}}*/
  /*io related{{{*/

  virtual std::string getVTKPointOrder() const {
    std::string content;
    for(size_t i = 0; i != getNumberOfNeighbors(0); ++i)
      content += std::to_string(getPoint(i).getIndex()) + " ";
    return content;
  }
  /*{{{*/
  /**     * Linear cells
   *   VTK_EMPTY_CELL       = 0,
   *   VTK_VERTEX           = 1,
   *   VTK_POLY_VERTEX      = 2,
   *   VTK_LINE             = 3,
   *   VTK_POLY_LINE        = 4,
   *   VTK_TRIANGLE         = 5,
   *   VTK_TRIANGLE_STRIP   = 6,
   *   VTK_POLYGON          = 7,
   *   VTK_PIXEL            = 8,
   *   VTK_QUAD             = 9,
   *   VTK_TETRA            = 10,
   *   VTK_VOXEL            = 11,
   *   VTK_HEXAHEDRON       = 12,
   *   VTK_WEDGE            = 13,
   *   VTK_PYRAMID          = 14,
   *   VTK_PENTAGONAL_PRISM = 15,
   *   VTK_HEXAGONAL_PRISM  = 16,
   *
   *   * Quadratic, isoparametric cells
   *   VTK_QUADRATIC_EDGE                   = 21,
   *   VTK_QUADRATIC_TRIANGLE               = 22,
   *   VTK_QUADRATIC_QUAD                   = 23,
   *   VTK_QUADRATIC_POLYGON                = 36,
   *   VTK_QUADRATIC_TETRA                  = 24,
   *   VTK_QUADRATIC_HEXAHEDRON             = 25,
   *   VTK_QUADRATIC_WEDGE                  = 26,
   *   VTK_QUADRATIC_PYRAMID                = 27,
   *   VTK_BIQUADRATIC_QUAD                 = 28,
   *   VTK_TRIQUADRATIC_HEXAHEDRON          = 29,
   *   VTK_QUADRATIC_LINEAR_QUAD            = 30,
   *   VTK_QUADRATIC_LINEAR_WEDGE           = 31,
   *   VTK_BIQUADRATIC_QUADRATIC_WEDGE      = 32,
   *   VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON = 33,
   *   VTK_BIQUADRATIC_TRIANGLE             = 34,
   *
   *   * Cubic, isoparametric cell
   *   VTK_CUBIC_LINE                       = 35,
   *
   *   * Special class of cells formed by convex group of points
   *   VTK_CONVEX_POINT_SET = 41,
   *
   *   * Polyhedron cell (consisting of polygonal faces)
   *   VTK_POLYHEDRON = 42,
   *
   *   * Higher order cells in parametric form
   *   VTK_PARAMETRIC_CURVE        = 51,
   *   VTK_PARAMETRIC_SURFACE      = 52,
   *   VTK_PARAMETRIC_TRI_SURFACE  = 53,
   *   VTK_PARAMETRIC_QUAD_SURFACE = 54,
   *   VTK_PARAMETRIC_TETRA_REGION = 55,
   *   VTK_PARAMETRIC_HEX_REGION   = 56,
   *
   *   * Higher order cells
   *   VTK_HIGHER_ORDER_EDGE        = 60,
   *   VTK_HIGHER_ORDER_TRIANGLE    = 61,
   *   VTK_HIGHER_ORDER_QUAD        = 62,
   *   VTK_HIGHER_ORDER_POLYGON     = 63,
   *   VTK_HIGHER_ORDER_TETRAHEDRON = 64,
   *   VTK_HIGHER_ORDER_WEDGE       = 65,
   *   VTK_HIGHER_ORDER_PYRAMID     = 66,
   *   VTK_HIGHER_ORDER_HEXAHEDRON  = 67,
   *
   *   * Arbitrary order Lagrange elements (formulated separated from generic
   *   higher order cells) VTK_LAGRANGE_CURVE           = 68,
   *   VTK_LAGRANGE_TRIANGLE        = 69,
   *   VTK_LAGRANGE_QUADRILATERAL   = 70,
   *   VTK_LAGRANGE_TETRAHEDRON     = 71,
   *   VTK_LAGRANGE_HEXAHEDRON      = 72,
   *   VTK_LAGRANGE_WEDGE           = 73,
   *   VTK_LAGRANGE_PYRAMID         = 74,
   */
  /*}}}*/
  virtual size_t getVTKTypeID() const {
    GEOXFEM_NOT_IMPLEMENTED();
    return 0;
  }
  virtual void
  plot(io::DataPloter&    plt,
       const std::string& seg_cfg = " w l ls 1",
       const std::string& gp_cfg  = " w p pt 0 lc rgb \"blue\"") const {
    for(auto& segment : segmentsIn(*this))
      segment->plot(plt, seg_cfg);
    if(!gp_cfg.empty())
      for(auto& gp : gpsIn(*this))
        gp.plot(plt, gp_cfg);
    return;
  }
  /**
   * @brief operator<<
   *
   * @param os
   * @param shp
   *
   * @return
   */
  friend std::ostream& operator<<(std::ostream& os, const Shape& shp) {
    os << shp.getClassName() << " #" << shp.getIndex() << '(';
    for(auto& vertex : neighborsIn(shp, 0))
      os << vertex->getIndex() << ",";
    os << ") with center: " << getCentroid(shp);
    return os;
  }
  /*}}}*/
  /*sub shape related{{{*/

  /**
   * @brief hasParent
   *
   * @return
   */
  bool hasParent() const {
    return m_parent != nullptr;
  }
  /**
   * @brief getParent
   *
   * @return
   */
  Shape* getParent() const {
    return m_parent;
  }
  /**
   * @brief setParent
   *
   * @param parent
   */
  void setParent(Shape* parent) {
    m_parent = parent;
  }
  /**
   * @brief Check if the object has a child.
   *
   * @return
   */
  bool hasChild() const {
    return !m_list_child.empty();
  }
  /**
   * @brief get number of child objects
   *
   * @return
   */
  size_t getNumberOfChild() const {
    return m_list_child.size();
  }
  /**
   * @brief clearChild
   */
  void clearChild() {
    for(auto v : m_list_child)
      delete v;
    m_list_child.clear();
  }
  /**
   * @brief Set child
   *
   * @param _pChild
   */
  void addChild(SubShape* child) {
    child->setParent(this);
    m_list_child.emplace_back(child);
  }
  /**
   * @brief deleteChild
   *
   * @param child
   */
  void deleteChild(SubShape* child) {
    auto it = std::find(m_list_child.begin(), m_list_child.end(), child);
    delete *it;
    m_list_child.erase(it);
    return;
  }
  /**
   * @brief getChildren
   *
   * @return
   */
  SubShapeList& getChildren() {
    return m_list_child;
  }
  /*}}}*/

private:
  GEOXFEM_Hash    m_hashkey;
  bool            m_boundary;
  GaussList       m_list_gauss;
  const Material* m_material;
  bool            m_xfem;
  Shape*          m_parent;
  SubShapeList    m_list_child;

protected:
  std::vector<ShapeList> m_listconnectivity;
};
}  // namespace geometry

namespace traits {

template <typename dimensionType, typename degree>
struct dimension_type</*{{{*/
                      geometry::ShapeInterface<dimensionType, degree>> {
  using type = dimensionType;
}; /*}}}*/

template <typename dimensionType, typename degree>
struct dimension</*{{{*/
                 geometry::ShapeInterface<dimensionType, degree>> {
  static constexpr size_t value = dimension<dimensionType>::value;
}; /*}}}*/

template <typename DimensionType, typename degree>
struct get_topology_dimension</*{{{*/
                              geometry::ShapeInterface<DimensionType, degree>> {
  using Shape = geometry::ShapeInterface<DimensionType, degree>;
  static size_t apply(const Shape& shp) {
    return shp.getTopologyDimension();
  }
}; /*}}}*/

}  // namespace traits

namespace algorithm {

/**
 * @brief
 */
template <typename dimensionType, typename degree>
struct clock_wise<geometry::ShapeInterface<dimensionType, degree>> { /*{{{*/
  using Shape      = geometry::ShapeInterface<dimensionType, degree>;
  using point_type = typename Shape::point_type;
  point_type m_centre;
  clock_wise(const point_type& _centre) : m_centre(_centre) {}
  template <typename Geometry>
  clock_wise(const Geometry& geometry) {
    m_centre = geometry.getCentroid();
  }
  bool operator()(const Shape& shp1, const Shape& shp2) const {
    point_type  p1  = shp1.getHeadPoint();
    point_type  p2  = shp2.getHeadPoint();
    const _real th1 = atan2(p1.y() - m_centre.y(), p1.x() - m_centre.x());
    const _real th2 = atan2(p2.y() - m_centre.y(), p2.x() - m_centre.x());
    return th1 > th2;
  }
  bool operator()(point_type* shp1, point_type* shp2) const {
    return operator()(*shp1, *shp2);
  }
}; /*}}}*/

/** get integration value through all gauss points in current element
 *
 * @param func
 * @return
 */
template <typename DimensionType, typename degree, typename Func>
auto integral(/*{{{*/
              const geometry::ShapeInterface<DimensionType, degree>& shp,
              const Func&                                            func)
  -> decltype(func(shp, model::GaussPoint<DimensionType>{})) {
  using Gauss        = model::GaussPoint<DimensionType>;
  Gauss          gp1 = shp.getGaussPoint(0);
  decltype(auto) ans = func(shp, gp1);
  ans *= 0;
  for(const auto& gp : gpsIn(shp))
    ans += func(shp, gp) * gp.getWeight();
  return ans;
} /*}}}*/
template <typename DimensionType, typename degree, typename Func>
auto integral(/*{{{*/
              const geometry::ShapeInterface<DimensionType, degree>& shp,
              const Func&&                                           func)
  -> decltype(func(shp, model::GaussPoint<DimensionType>{})) {
  using Gauss        = model::GaussPoint<DimensionType>;
  Gauss          gp1 = shp.getGaussPoint(0);
  decltype(auto) ans = func(shp, gp1);
  ans *= 0;
  for(const auto& gp : gpsIn(shp))
    ans += func(shp, gp) * gp.getWeight();
  return ans;
} /*}}}*/
/** get weighted average value through all gauss points in current element
 *
 * @param func
 * @return
 */
template <typename DimensionType, typename degree, typename Func>
auto weighted_average(/*{{{*/
                      const geometry::ShapeInterface<DimensionType, degree>&
                                   shp,
                      const Func&& func)
  -> decltype(func(shp, model::GaussPoint<DimensionType>{})) {
  return integral(shp, func) / shp.getArea();
} /*}}}*/
template <typename DimensionType, typename degree, typename Func>
auto weighted_average(/*{{{*/
                      const geometry::ShapeInterface<DimensionType, degree>&
                                  shp,
                      const Func& func)
  -> decltype(func(shp, model::GaussPoint<DimensionType>{})) {
  return integral(shp, func) / shp.getArea();
} /*}}}*/
/** get value throught interpolation
 *
 * @param gcoord
 * @param func    :
 * @return
 */
template <typename Shape, typename VertexFunc>
auto interpolate(/*{{{*/
                 const Shape&      shp,
                 const VectorXR&   gcoord,
                 const VertexFunc& func) -> decltype(auto) {
  decltype(auto) ans = func(*shp.getNeighbourBegin(0)) * 0;
  VectorXR       N   = shp.getStdNByGlobal(gcoord);
  size_t         i   = 0;
  for(auto& v : neighborsIn(shp, 0))
    ans += func(&(v->toVertex())) * N[i++];
  return ans;
} /*}}}*/

/**
 * @brief G2V_mapping
 *
 * @param shp
 */
template <typename DimensionType, typename degree>
void G2V_mapping(/*{{{*/
                 const geometry::ShapeInterface<DimensionType, degree>& shp) {
  for(auto& gp : gpsIn(shp)) {
    VectorXR N      = shp.getStdNByGlobal(gp.getLocalCoord());
    size_t   loc    = 0;
    _real    weight = N[loc++] * gp.getWeight();
    for(auto& vertex : neighborsIn(shp, 0)) {
      auto& v_gp     = vertex->toVertex().getGauss();
      _real v_weight = v_gp.getWeight();
      v_gp += gp * weight / v_weight;
    }
  }
  return;
} /*}}}*/

}  // namespace algorithm

namespace model {

template <typename CasePolicy>
struct default_gauss_num { /*{{{*/
  using ptree = boost::property_tree::ptree;
  static std::unordered_map<std::string, size_t>&
  set_get(const ptree& data = ptree{}) {
    static std::unordered_map<std::string, size_t> pt = {
      {"segment", 3}, {"triangle", 3},  {"quadrangle", 4},
      {"tetra", 4},   {"tri_prism", 6}, {"octahedron", 8}};
    if(!data.empty())
      for(auto& gp_num_data : data) {
        std::string name   = gp_num_data.second.get<std::string>("name");
        size_t      num_gp = gp_num_data.second.get<size_t>("value");
        pt[CasePolicy::apply(name)] = num_gp;
      }
    return pt;
  }
};  // namespace model
using i_default_gauss_num = default_gauss_num<util::ignore_case>; /*}}}*/

template <typename DimensionType, typename degree>
arithmetic::Stress<DimensionType>
GetCurrSigma(/*{{{*/
             const geometry::ShapeInterface<DimensionType, degree>& shp) {
  using Shape = geometry::ShapeInterface<DimensionType, degree>;
  using Gauss = model::GaussPoint<DimensionType>;

  return algorithm::weighted_average(
    shp, [](const Shape&, const Gauss& gp) { return gp.getCurrSigma(); });
} /*}}}*/
template <typename DimensionType, typename degree>
arithmetic::Strain<DimensionType>
GetCurrEpsilon(/*{{{*/
               const geometry::ShapeInterface<DimensionType, degree>& shp) {
  using Shape = geometry::ShapeInterface<DimensionType, degree>;
  using Gauss = model::GaussPoint<DimensionType>;

  return algorithm::weighted_average(
    shp, [](const Shape&, const Gauss& gp) { return gp.getEpsilon(); });
} /*}}}*/

template <typename DimensionType, typename degree>
void add_gauss(GaussPoint<DimensionType>&                       gp,
               geometry::ShapeInterface<DimensionType, degree>& shp) {
  gp = algorithm::interpolate(
    shp, gp.getCoord(), [](const auto* v) { return v->toVertex().getGauss(); });
  shp.append(gp);
  return;
}

template <typename DimensionType, typename degree>
void buildGaussPoints(/*{{{*/
                      geometry::ShapeInterface<DimensionType, degree>& shp) {
  shp.clearGaussPoints();
  std::vector<GaussPoint<DimensionType>> gp_list;
  if(shp.hasChild())
    for(auto child : childrenIn(shp))
      child->getStdGaussPoints(gp_list);
  else
    shp.getStdGaussPoints(gp_list);
  for(GaussPoint<DimensionType>& gp : gp_list)
    add_gauss(gp, shp);
  return;
} /*}}}*/

}  // namespace model

namespace io {

struct vtu_cell_offset { /*{{{*/
  template <typename DimensionType, typename degree>
  std::string
  operator()(const geometry::ShapeInterface<DimensionType, degree>& shp) const {
    static size_t np = 0;
    if(shp.getIndex() == -1) {
      np *= 0;
      return "\0";
    }
    np += shp.getNumberOfPoint();
    return std::to_string(np);
  }
};                             /*}}}*/
struct vtu_cell_connectivity { /*{{{*/
  template <typename DimensionType, typename degree>
  std::string
  operator()(const geometry::ShapeInterface<DimensionType, degree>& shp) const {
    std::string content;
    for(auto& v : neighborsIn(shp, 0))
      content += std::to_string(v->getIndex()) + " ";
    return content;
  }
};                     /*}}}*/
struct vtu_cell_type { /*{{{*/
  template <typename DimensionType, typename degree>
  std::string
  operator()(const geometry::ShapeInterface<DimensionType, degree>& shp) const {
    return std::to_string(shp.getVTKTypeID());
  }
};                      /*}}}*/
struct vtu_cell_sigma { /*{{{*/
  template <typename DimensionType, typename degree>
  std::string
  operator()(const geometry::ShapeInterface<DimensionType, degree>& shp) const {
    auto stress = model::GetCurrSigma(shp);
    /* Create an output string stream */
    std::ostringstream streamObj;
    /* Add double to stream */
    streamObj << stress(0, 0) << " "   //
              << stress(1, 0) << " "   //
              << stress(2, 0) << "\n"  //
              << stress(0, 1) << " "   //
              << stress(1, 1) << " "   //
              << stress(2, 1) << "\n"  //
              << stress(0, 2) << " "   //
              << stress(1, 2) << " "   //
              << stress(2, 2) << "\n";
    /* Get string from output string stream */
    return streamObj.str();
  }
};                        /*}}}*/
struct vtu_cell_epsilon { /*{{{*/
  template <typename DimensionType, typename degree>
  std::string
  operator()(const geometry::ShapeInterface<DimensionType, degree>& shp) const {
    auto strain = model::GetCurrEpsilon(shp);
    /* Create an output string stream */
    std::ostringstream streamObj;
    /* Add double to stream */
    streamObj << strain(0, 0) << " "   //
              << strain(1, 0) << " "   //
              << strain(2, 0) << "\n"  //
              << strain(0, 1) << " "   //
              << strain(1, 1) << " "   //
              << strain(2, 1) << "\n"  //
              << strain(0, 2) << " "   //
              << strain(1, 2) << " "   //
              << strain(2, 2) << "\n";
    /* Get string from output string stream */
    return streamObj.str();
  }
}; /*}}}*/
}  // namespace io

}  // namespace geoxfem

struct shapeHasher { /*{{{*/
  template <typename DimensionType, typename degree>
  size_t operator()(const std::unique_ptr<
                    geoxfem::geometry::ShapeInterface<DimensionType, degree>>&
                      _spShape) const {
    return std::hash<std::string>{}(_spShape->getHashKey());
  }
  template <typename DimensionType, typename degree>
  size_t operator()(
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>* p_shp)
    const {
    return std::hash<std::string>{}(p_shp->getHashKey());
  }
  template <typename DimensionType, typename degree>
  size_t operator()(
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>& shp) const {
    return std::hash<std::string>{}(shp.getHashKey());
  }
};                /*}}}*/
struct shapeCmp { /*{{{*/
  template <typename DimensionType, typename degree>
  bool operator()(
    const std::unique_ptr<
      geoxfem::geometry::ShapeInterface<DimensionType, degree>>& lhs,
    const std::unique_ptr<
      geoxfem::geometry::ShapeInterface<DimensionType, degree>>& rhs) const {
    return lhs->getHashKey() == rhs->getHashKey();
  }
  template <typename DimensionType, typename degree>
  bool operator()(
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>* lhs,
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>* rhs) const {
    return lhs->getHashKey() == rhs->getHashKey();
  }
  template <typename DimensionType, typename degree>
  bool operator()(
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>& lhs,
    const geoxfem::geometry::ShapeInterface<DimensionType, degree>& rhs) const {
    return lhs.getHashKey() == rhs.getHashKey();
  }
}; /*}}}*/

#endif /* I_SHAPE_H */
